import matplotlib
#matplotlib.use('WXAgg')
from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *

#from brian.globalprefs import *
#set_global_preferences(useweave=True,usecodegen=True,usecodegenweave=True,gcc_options=['-ffast-math','-march=native'],usenewpropagate=True)
import mydataMN as mydatmn
md = mydatmn.MNVars()
#print md.DUR

defaultclock.dt=0.0025*ms


def i_inj(t):
    if t%(10*ms) <= 1*ms and t%(10*ms)>= 0 *ms:
        return 0.5 *nA
    else:
        return 0*nA
    

locDUR =500


def spToList(a):
    l = list()
    
    for x in a.spiketimes:
        l+=((a[x]/20*1000/md.BINSIZE)/ms).tolist()
    return l

def napvtrap(x,y):
    a=x/y
    if abs(a)<0.0001:
        return  1/(2+a)
    else:
        return 1/(exp(a)+1)
    #return 1/(exp(a)+1)
    
def napefun(x,y):
    a = x/y 
    if abs(a)<0.0001:
        return 1+a
    else:
        return exp(a)
    #return exp(a)
            
 
def m_inf_na(V):
    return napvtrap(-(V+35),7.8)   
    
def h_inf_na(V):
    return napvtrap((V+55),7)

def tau_h_na(V):
    return msecond*(30/(napefun((V+50),15)+napefun(-(V+50),16)))

def m_inf_k(V):
    return napvtrap(-(V+28),15)

def tau_m_k(V):
    return msecond*(7/(napefun(V+40,40)+napefun(-(V+40),50)))

def m_inf_can(V):
    return napvtrap(-(V+30),5)

def h_inf_can(V):
    return napvtrap(V+45,5)

def m_inf_cal(V):
    return napvtrap(-(V+40),7)


def m_inf_nap(V):
    return napvtrap(-V-47.1,3.1)

def h_inf_nap(V):
    return napvtrap(V+59,8)

def tau_h_nap(V):
    return msecond*(2*1200/(napefun(V+59,16)+napefun(-V-59,16)))

def current_inj(t):
    if t > 0*msecond:
        return -0.5*namp
    else:
        return 0*mamp

# The model

#eqsPF+=OrnsteinUhlenbeck('I',mu=0.005*nA,sigma=0.002*nA,tau=10*ms)

# soma: INa,IK,ICaN,IKCa,IL,IC
# dendrite:
eqsMN=Equations('''
dvm/dt = (i_inj(t)-(md.gna*mna*mna*mna*hna*(vm-md.Ena))-(md.gk*mk*mk*mk*mk*(vm-md.Ek))-(md.gcan*mcan*mcan*hcan*(vm-md.eCa))-(md.gkca*ca/(ca+md.kd)*(vm-md.Ek))-(md.gl*(vm-md.Elmn))-(md.gc/md.p*(vm-vmd)))/md.cpf : mvolt
dmna/dt = (m_inf_na(vm/mvolt)-mna)/md.taumna : 1
dhna/dt = (h_inf_na(vm/mvolt)-hna)/tau_h_na(vm/mvolt) : 1
dmk/dt = (m_inf_k(vm/mvolt)-mk)/tau_m_k(vm/mvolt) : 1
dmcan/dt = (m_inf_can(vm/mvolt)-mcan)/md.tau_m_can : 1
dhcan/dt = (h_inf_can(vm/mvolt)-hcan)/md.tau_h_can : 1
dca/dt = md.f*(-md.alpha*(md.gcan*mcan*mcan*hcan*(vm-md.eCa))-md.kca*ca) : mole*dmetre**-3

dvmd/dt = (-(md.gnap*mnapd*hnapd*(vmd-md.Ena))-(md.gcan*mcand*mcand*hcand*(vmd-md.eCa))-(md.gcal*mcald*(vmd-md.eCa))-(md.gkca*cad/(cad+md.kd)*(vmd-md.Ek))-(md.gl*(vmd-md.Elmn))-(md.gc/(1-md.p)*(vmd-vm)))/md.cpf :mvolt
dmnapd/dt = (m_inf_nap(vmd/mvolt)-mnapd)/md.tau_m_nap : 1
dhnapd/dt = (h_inf_nap(vmd/mvolt)-hnapd)/tau_h_nap(vmd/mvolt) : 1
dmcand/dt = (m_inf_can(vmd/mvolt)-mcand)/md.tau_m_can : 1
dhcand/dt = (h_inf_can(vmd/mvolt)-hcand)/md.tau_h_can : 1
dmcald/dt = (m_inf_cal(vmd/mvolt)-mcald)/md.tau_m_cal : 1
dcad/dt = md.f*(-md.alpha*((md.gcan*mcand*mcand*hcand*(vmd-md.eCa))+(md.gcal*mcald*(vmd-md.eCa)))-md.kca*cad) : mole*dmetre**-3


''')
#dvmd/dt = (-(md.gnap*mnap*hnap*(vmd-md.Ena))-(md.gcan*mcan*mcan*hcan*(vmd-md.eCa)-(gkca*ca/(ca+0.2)*(vmd-eK)-md.glmn*(vmd-md.El))
eqsMN+=OrnsteinUhlenbeck('I',mu=0.0*nA,sigma=0.025*nA,tau=10*ms)

#FREQ=30*Hz
#spiketimes = [(0,5*ms),(1,5*ms),(2,5*ms),(3,5*ms),(4,5*ms),(5,5*ms),(6,5*ms),(7,5*ms),(8,5*ms),(9,5*ms),(10,5*ms),(11,5*ms),(12,5*ms),(13,5*ms),(14,5*ms),(15,5*ms),(16,5*ms),(17,5*ms),(18,5*ms),(19,5*ms),]
print eqsMN
    
#SG = SpikeGeneratorGroup(20,spiketimes,period=1/FREQ)

MN=NeuronGroup(1,model=eqsMN,
threshold=EmpiricalThreshold(threshold=-20*mV),
implicit=True,method='exponential_Euler')


# Initialization



vin =-61.4557
MN.vm=vin*mV
MN.mna = m_inf_na(vin)
MN.hna = h_inf_na(vin)
MN.mk = m_inf_k(vin)
MN.mcan = m_inf_can(vin)
MN.hcan = h_inf_can(vin)
MN.ca = -md.alpha*(md.gcan*MN.mcan*MN.mcan*MN.hcan*(MN.vm-md.eCa))/md.kca

vin = -59.835
MN.vmd=vin*mV
MN.mnapd = m_inf_nap(vin)
MN.hnapd = h_inf_nap(vin)
MN.mcand = m_inf_can(vin)
MN.hcand = h_inf_can(vin)
MN.mcald = m_inf_cal(vin)
MN.cad = -md.alpha*((md.gcan*MN.mcand*MN.mcand*MN.hcand*(MN.vmd-md.eCa))+(md.gcal*MN.mcald*(MN.vmd-md.eCa)))/md.kca


#IN.El = -57.5*mV * randn(len(P))*2.875*mV

# Record a few trace
tvm=StateMonitor(MN,'vm',record=[0])
tvmd=StateMonitor(MN,'vmd',record=[0])

tmna=StateMonitor(MN,'mna',record=[0])
thna=StateMonitor(MN,'hna',record=[0])
tmk=StateMonitor(MN,'mk',record=[0])
tmcan=StateMonitor(MN,'mcan',record=[0])
thcan=StateMonitor(MN,'hcan',record=[0])
tca=StateMonitor(MN,'ca',record=[0])

#spikesIN= SpikeMonitor(IN)

run(locDUR*msecond)
figure(1)
#plot(tracePFgi.times/ms,tracePFgi[0]/nA)
#plot(tracePFge.times/ms,tracePFge[0]/nA)

plot(tvmd.times/ms,tvmd[0]/mV)

#subplot(711)
#plot(tvm.times/ms,tvm[0]/mV)
#subplot(712)
#plot(tmna.times/ms,tmna[0])
#subplot(713)
#plot(thna.times/ms,thna[0])
#subplot(714)
#plot(tmk.times/ms,tmk[0])
#subplot(715)
#plot(tmcan.times/ms,tmcan[0])
#subplot(716)
#plot(thcan.times/ms,thcan[0])
#subplot(717)
#plot(tca.times/ms,tca[0]/mm)

figure(2)
plot(tvm.times/ms,tvm[0]/mV)

#subplot(412)
#x = spToList(spikesIN)
#hist(x,floor(md.DUR/md.BINSIZE),range=[0,md.DUR],histtype='step')
show()